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FOREWORD 

Ttvjg  report  was  prepared  by  the  Directorate  of  Chemical  Sciences, 
Frank  J.  Seiler  Research  Laboratory,  United  States  Air  Force  Academy, 
Colorado.  The  work  has  initiated  under  Project  No.  2303,  "Chemistry," 
Task  NO.  2303-F2,  "Physical  Chemistry  and  Electrochemistry,"  Work  Unit 
No.  2303-F2-06,  "Physical  and  Electrochemical  Measurements." 

Hie  report  covers  work  conducted  from  April  1975  to  December  1976. 
The  nanuscript  was  released  by  the  authors  for  publication  in  December 
1976. 

This  Technical  Report  has  been  reviewed  and  approved. 
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INTRODUCTION 


The  customary  procedure  followed  in  least  squares  fitting  of 
experimental  data  to  an  assumed  equation  is  to  consider  all  of  the 
errors  to  be  in  one  of  the  variables,  and  that  all  other  variables  are 
exact.  For  example,  consider  a least  squares  fit  of  the  equation 
y = a + bx 

Typically,  x is  considered  to  be  an  independent  variable,  not  subject  to 
error,  and  y is  considered  to  be  a dependent  variable.  If  yQ  are  the 
observed  values  of  the  dependent  variable,  we  seek  to  minimize  with 
respect  to  a and  b the  quantity  Z (yQ-y)  • 

While  one  may  indeed  mathematically  consider  x an  independent  variable 
and  not  subject  to  error,  and  y to  contain  the  error,  one  might  equally 
well  reverse  the  roles  of  x and  y.  In  fact,  read  physical  situations 
usually  involve  error  in  both  x and  y,  and  it  may  not  always  be  stated 
that  one  variable  is  more  suitable  for  use  as  an  independent  variable 
than  the  other. 

It  is  possible  to  treat  both  x and  y as  subject  to  error  in  least 
squares  fitting,  but  it  does  not  often  seem  to  be  done.  We  describe 
here  a technique  for  least  squares  fitting  which  minimizes  the  distance 
from  the  experimental  point  to  the  assured  equation.  The  technique  is 
applicable  to  any  number  of  variables  subject  to  error  and  any  number  of 
adjustable  constants  in  the  assumed  equation. 


DISCUSSION 


Consider  the  curve 
y = f (x) 

in  Figure  1.  The  minimum  distance  from  the  experimental  point  (x,y)  to 
the  assumed  equation  is  distance  d in  Figure  1 (a) . Providing  the 
curvature  is  not  too  great,  or  expressed  alternately,  providing  Ax  and  Ay 
are  sufficiently  snail,  the  distance  d may  be  approximated  by  the 
perpendicular  distance  z from  point  (x,y)  to  the  chord  (x^1’)  (x',y)  in 
Figure  1 (b) . Obviously  z may  be  snaller  or  larger  than  d,  depending  on 
the  sign  of  the  curvature,  or  on  vrfiich  side  of  the  equation  point  (x,y) 
lies. 

The  least  squares  technique  described  here  minimizes  the  quantity 

Ez2. 

The  discussion  below  will  treat  as  an  example  the  two-dimensional 
rapo  of  Figure  1(b)  with  a set  of  parallel  equations  for  the  n-dimensional 
case.  Chord  (x,y*)  (x',y)  may  be  considered  as  a one-dimensional  surface. 
The  general  equation  for  a surface  is 
Ax  + By  + C = 0 
n 

E a.x.  = 0 where  x = 1 
i=0  1 1 

Let  (x,y) [(x1,x2,x3,...xn)]  be  an  experimental  point,  and  x'  and  y' 

fx,  ' ,x~',x  ',...x  ']  be  the  values  of  the  variables  calculated  from  the 
1 z 3 n 

assumed  equation  and  the  appropriate  remaining  independent  variables  for 
the  given  experimental  point.  That  is 
x'  = f(y  and  the  parameters) 
y'  = f (x  and  the  parameters) 


T"  + A“  + C = 0 
Ax  Ay 

c = -(£-+*-) 
u ’•Ax  Ay-1 

n 

a + a.x.'  + l a.x.  = 0 

° 3 3 irtfirt  1 1 

ao  * -ajxj'  I.  ai*i 


From  analytic  geometry 

2 _ (Ax  + By  + C)' 

2 2 "2 
A2  + B 


2 la°  +iiaiXir 

Z = ~ 2 

I <ai>2 

i=l 


Substituting  Equations  (7)  and  (8)  into  Equation  (9)  , we  have 


if  (i)2  + 


VM’2 

Unlike  z^,  the  distance  z will  have  a sign,  which  must  be  selected 
when  the  square  root  is  taken.  Since  z is  a vector  quantity,  its  magnitude, 
which  corresponds  to  the  positive  root,  is  used  in  the  fitting  procedure. 
The  procedure  produces  a least  squares  sun  of  these  magnitudes.  Signs 
for  the  corrections  to  the  parameters  arise  from  the  partial  derivatives 
of  z with  respect  to  the  parameters  and  not  from  the  sign  of  z itself. 

The  assumed  functional  relationship  among  the  parameters  will  in 
general  involve  one  or  more  parameters,  or  adjustable  constants,  the 


5 
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values  of  which  we  must  find  to  minimize  the  value  Xz  . Since  we 
the  true  equation  is  very  closely  approximated  by  the  surface  of 
Equation  (3) , we  may  write  for  each  of  the  parameters 


9z  _3z  , 3x'(V,  . 3z  , 


30k”  3Ax  l" 


36, 


■J  + w h 


36, 


-) 


n . 3x. ' (0,  ) 

oZ  r ( oZ  r 1 iC  n A 

30k  l3Ax7  1 m. 


36, 


Equation  (11)  sinplifies  to 


3z  _ _3  r 1 ax'(ek)  , 1 

"k  (Ax)3  "k 


3y'(ek) 

36,. 


(11) 

dl') 


(12) 


3z  _ _3  " r 1 3xi'  (ek\ 

gg  — — z L *■  Y * 


36, 


(12') 


Jk  i=l  (Ax±)  k 

The  function  to  be  minimized,  given  by  Equation  (10) , and  the 


derivatives  with  respect  to  each  of  the  parameters,  given  by  Equation  (12) , 
must  be  calculated  for  each  experimental  point,  and  are  the  quantities 
actually  used  in  the  least  squares  fitting  routines.  The  latter  are 


given  in  the  next  section  "Mathematical  Procedure." 

It  is  better  not  to  use  the  set  of  variables,  {x^ , in  the  actual 
units  with  which  they  were  measured,  for  the  choice  of  units  will  affect 
the  final  values  of  the  parameters,  {0k}.  Dimensionless  variables  - 
independent  of  the  choice  of  units  - may  be  obtained  by  introducing  the 
variables  in  terms  of  the  estimated  uncertainty  or  experimental  error  in 
that  variable.  For  example,  if  a particular  variable,  Kt  may  be 
determined  within  an  uncertainty  of  6£,  the  dimensionless  variable, 

£*  = £/<$£#  is  used  in  Equations  (3)  through  (12).  Equations  (10)  and 


(12)  become,  respectively 


6 


Note  that  the  fitted  values  of  {0^}  will  still  be  in  «any>  units  as  would 
have  been  obtained  fran  using  Equations  (10)  and  (12) . 

Inspection  of  Equations  (13)  and  (14)  reveals  that  if  all  6xi  but 
one  are  made  extremely  snail,  the  minimum  distance  fit  reduces  to  a 
conventional  least  squares  fit  with  all  of  the  error  cast  into  the 
remaining  variable.  A useful  application  of  this  property  would  be 
when  the  experimenter  was  not  sure  of  the  data  measurement  errors  in 
one  or  more  of  the  variables,  and  wished  to  selectively  include  or 
exclude  individual  variables  from  fits. 


The  quantities  {Ax^}  = {x^  - xV } were  individually  calculated  by 
assuming  all  of  the  error  was  cast  into  that  particular  variable. 

Clearly  the  actual  error  in  the  variables  {x^}  is  less  than  {Ax^} . A nei 
set  of  quantities,  {A'x^} , is  defined  with  the  aid  of  Figure  2(a).  The 
assumption  is  again  made  that  the  true  locus  of  the  equation  lies 
negligibly  far  away  from  the  surface,  i.e. , d=z. 


From  analytic  geometry 
_2 


A'x.  = -r- 
i Ax. 
i 


A'x.*  = t— j 

l Ax^ 


Where  A'x^*  = A'x^/fix^.  If  the  original  error  estimate,  <5x^,  was 
correct,  then  ideally  A'x.  and  fix.  should  be  nearly  the  same  size. 


After  the  point  (xg,yg)  [xlg,  x2g,  x3s' • • has  been  identified 
it  is  possible  to  determine  the  validity  of  the  assurptions  that  the 
surface  and  the  true  locus  were  very  close.  It  may  be  seen  in  Figure  2(b) 
that  d-z  = z' , and  z'  nay  be  obtained  by  application  of  Equation  (10)  to 
the  point  (xs,y8)  [xlg,  x^. 


MATHEMATICAL  PROCEDURE 


Physical  measurements  generally  involve  the  simultaneous  measure- 
ment of  a nuifcer  (n)  of  related  physical  properties.  If  the  values  of 
these  properties  are  arranged  in  ordered  sets  of  n values,  vectors  or 
state  points  {xn}  are  obtained.  If  a functional  relationship  exists 
among  these  n variables,  containing  in  addition  to  the  variables  a 
number  (p)  of  constants  characteristic  of  the  physical  system,  then  a 
functional  surface,  f,  may  be  defined  such  that  f is  identically  zero. 

0 = f({xn^V) 

This  ordered  set  of  constants  |0p},  called  parameters,  which 
characterizes  the  system  also  forms  a vector.  For  instance,  given  the 
typical  vapor  pressure  equation,  0 = log  (P)  - + B) , the  ordered  set  of 

constants  (A,B)  are  sufficient  to  define  the  pressure,  P,  at  any 
temperature , T,  or  the  temperature  at  any  pressure.  However,  each 
different  physical  system  obeying  the  same  relationship  will  have  a 
different  set  of  parameters  (A,B) . 

In  principle  the  evaluation  of  these  parameters  for  a system  should 
be  simple.  If  there  are  p parameters,  then  p state  points  are  measured 
producing  p equations  in  as  many  unknowns  which  may  be  solved  for  the 
parameters.  TWo  problems  immediately  present  themselves  when  an  attempt 
is  made  to  utilize  this  method.  First,  if  the  functioned  relationship 
is  not  linear  in  the  parameters  {0^} , then  solution  of  the  equations  may 
be  ambiguous,  difficult,  or  impossible.  Second,  since  no  measurement  is 
exact  but  generally  yields  only  a few  significant  figures,  the  function 
f will  not  yield  zero,  but  seme  small  non-zero  value,  d («z) , when  the 
values  of  any  measured  state  point  are  inserted  in  the  function. 


The  initial  difficulty  of  non-linearity  may  be  overcame  by 
linearizing  the  equations  for  { 0p} . (To  avoid  confusion  later,  vectors 
and  matrices  will  be  indicated  henceforth  by  the  use  of  boldface  instead 
of  the  subscript-in-braces  notation  used  earlier.)  Ibis  is  done  by 
expanding  the  equations  for  0 in  a Taylor  series  and  retaining  terms  in 
A0^  to  only  first  order.  Thus  for  each  state  vector  an  equation  exists 
of  the  form 

("A6i)+  ("a62^  + ’ + 30  ^-A0p)  = Z0-Z0+A0 

Since  zQ  = f = 0,  the  right  hand  side  of  Equation  (18)  reduces  to 

Z0  “ Z0+A0  = "Z0+A0 

When  reduced  to  matrix  form  these  equations  became 


[22  22  ...22) 

'■30.  30-  30  J 

12  p 


-A0, 


-A0, 


,-A0  , 

L PJ 


= ^~Z0+A0^ 


or  defining  {||-}  = P'  and  {— A0  } = A 


P'A  = 2 


0+A0 


Using  such  an  expansion  greatly  simplifies  the  solution  of  the 
simultaneous  equations,  but  introduces  limitations  on  the  applicability 
in  the  following  manner.  Since  0 is  solved  for  by  estimating  errors  in 
0 using  a linear  function,  the  solution  depends  on  the  initial  estimate 
of  0 and  the  curvature  of  the  error  surface  at  this  point.  If  the 
function  is  linear  in  0 there  is  no  curvature  and  any  initial  0 will 
converge  to  the  correct  solution.  In  regions  where  the  function  is 
linear  or  nearly  linear,  then  0 can  be  solved  for  in  one  or  more 


iterations.  Where  the  function  deviates  significantly  from  linearity 
in  0 the  solutions  may  occur  in  regions  where  they  are  physically 
unrealistic  (convergence  to  an  incorrect  solution)  or  the  i I’naar  approx- 
imation may  be  sufficiently  in  error  to  move  further  and  further  fran 
any  solution  with  each  iteration  (divergence) . For  a well  behaved  and 
at  least  piece-wise  continuous  function  then  a solution  generally  exists. 
For  such  functions  the  degree  of  deviation  from  linearity  in  0 in  the 
region  of  the  solution  will  determine  how  close  to  the  solution  the 
initial  estimate  of  0 must  be  for  convergence  to  the  correct  solution. 

The  second  difficulty,  that  associated  with  the  error  of  measure- 
ment is  treated  in  the  following  nanner. 

Measurements  are  made  of  many  more  than  p state  vectors.  Since,  in 
general,  no  two  state  vectors  will  give  the  same  value  for  z,  this 
provides  a set  of  k equations  (k>p)  which  are  inconsistent.  The  "degree" 
of  inconsistency  however  is  small,  that  is  there  will  be  a distribution 
of  z's  about  zero  which  is  assumed  to  be  Gaussian  or  normal  and  is  a 
function  of  the  errors  in  the  measurements  of  the  state  variables.  The 
set  of  "overdetermined"  parameters  (i.e.,  those  defined  by  all  k equations 
instead  of  only  the  minimum  subset  p of  them)  cure  solved  for  by  establish- 
ing a criterion  for  the  "best"  parameters.  Logically  the  best  set  of 
paraneters  would  be  those  predicted  to  be  most  probable  from  the  Gaussian 
distribution  of  the  errors,  z.  This  would  be  difficult  to  do  using 
functions  of  the  distribution  itself  and  such  a procedure  turns  out  to 
be  unnecessary.  It  has  been  shown  that  the  "least  squares  requirement" 
gives  the  most  probable  solution  for  the  parameters  (1) . This  "least 
squares  requirement"  therefore  is  used  ocnnonly  as  the  criterion  for  the 


: 


best  fit.  The  least  squares  requirement  is  very  sinple  in  concept.  If 
the  error  is  defined  in  terms  of  the  distance  of  the  state  vectors  from 
the  functional  surface,  the  error  surface,  T,  is  a function  of  0 defined 
by  the  sun  of  the  squares  of  all  individual  errors,  z.  This  sun  of  the 
squared  errors,  Ez  , is  minimized  by  adjustment  of  the  parameters  0. 

The  adjustment  is  accomplished  by  using  standard  techniques  from 
the  calculus  and  linear  algebra.  The  unconstrained  extremun  for  any 
well-behaved  function  which  is  linear  in  several  variables,  is  found  by 
setting  the  partial  derivatives  of  the  function  with  respect  to  each  of 
the  variables  equal  to  zero  and  solving  the  system  of  equations  obtained. 

These  cure  called  the  normal  equations.  The  error  surface  is  defined  to  be 

T = Ez2  (21) 

The  set  of  equations  obtained  by  setting  the  partial  derivatives  of  T 
with  respect  to  each  of  the  parameters  equal  to  zero  are  the  normal 
equations  for  the  system.  That  is 


It  can  be  shown  in  the  individual  cases  of  practical  interest  that 
the  extremum  so  obtained  is  a minimum  This  can  be  seen  from  the  fact 
that  for  any  functional  surface  a maximum  distance  from  the  surface  to 
the  fixed  set  of  data  points  does  not  exist,  i.e. , the  surface  can  always 
be  moved  further  away  from  all  the  points. 

Margulis  (2)  also  points  out  that  the  "least  squares  criterion" 
satisfies  the  following  intuitively  desirable  requirements  for  a "best  fit": 
(1)  the  error  is  expressed  as  an  unsigned  value  including  ALL 
measured  vector  errors. 
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(2)  the  overall  error  is  zero  only  if  ALL  errors  are  zero, 

(3)  the  overall  error  increases  if  any  ONE  error  increases. 

It  is  necessary  to  build  the  system  of  normal  equations  and  solve 

them.  The  system  could  be  constructed  by  evaluation  of  the  appropriate 
derivatives  of  T,  but  a sinpler  method  exists  which  has  been  shown  to  be 
equivalent.  This  latter  method  is  based  on  the  construction  of  the 
correction  Equation  (18) . For  each  measured  state  vector  x^ , one 
equation  like  Equation  (20)  is  obtained. 

Pj'A  = z.  (23) 

Hie  vector  A has  no  subscript  since  a ccnmon,  unique  set  of 
corrections  for  all  the  state  vectors  is  assumed  to  exist.  Since  each 
vector  Pj ' is  a row  of  ordered  values,  the  set  of  equations  for  k state 
vectors  will  form  a matrix  with  k rows  (one  for  each  state  vector)  ani 
p columns  (one  for  each  parameter) . 

The  column  vector  A is  to  be  evaluated  and  added  to  n0  to  obtain 
an  improved  0 . The  column  vector  Z formed  by  the  k values  of  the 


1 


Zj  is  used  in  the  solution.  When  all  of  the  elements  are  collected, 
the  set  of  correction  equations  appear  as: 
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or 

PA  = Z 

m 

Prenrultiplying  this  equation  on  both  sides  by  P the  transpose  of 
matrix  P gives  a matrix  equation  that  is  identical  to  the  normal 
equation  system. 

ptpa  =ptz 

If  the  products  P T/>  and  PTZ  are  assigned  the  symbols  H and  B 
respectively  then  the  matrix  equation 

HA  = B 

is  the  one  that  must  be  solved.  Inversion  of  the  square  matrix  H to 
give  //_1  provides  this  solution.  Premultiplying  both  sides  of 
Equation  (27)  by  H ^ obtains 

H-1H  A = H_1B 
by  definition  H ^ HA  - EA  - A 

where  E is  the  unit  matrix  and  therefore 

A = H-1B 

It  would  appear  at  this  point  that  the  simplification  introduced 
by  using  this  method  is  negligible,  but  the  construction  of  the  matrices 
H and  B is  particularly  simple. 

Initially  a null  natrix  (all  elements  zero)  for  both  H and  B is 
constructed.  //  is  a p x p matrix  (i.e. , p rows  and  p columns)  and 
B is  a p x 1 natrix. 

_» 

For  each  state  point  for  1 to  k the  vector  P described  earlier 
is  evaluated.  The  elements  of  this  vector  are  the  partials  of  z with 
respect  to  each  of  the  parameters  0^  (i  = 1 to  p) . z is  also  evaluated. 


The  elements  of  H are  then  augmented  by  adding  to  each  element 
hij  ^ product  pA’  x p.'.  Similarly  the  product  z x p.  is  added  to  b.. 

This  process  is  identical  to  the  following  matrix  operation. 

T 

^(augmented)  ” ^(unaugmented)  + ^ P ) ( P') 
and 

^(augmented)  “ ^(unaugmented)  + ^ ^ 

Once  these  matrices  have  been  constructed,  the  inversion  of  H may 
be  carried  out  using  any  standard  procedure. 

The  mathematical  discussion  in  this  section  has  been  phrased  in 
terms  of  the  variables  fx^}  and  the  minimum  distance,  z;  however,  the 
procedure  is  equally  valid  for  the  use  of  reduced  variables  {x>}  and  z*. 


APPLICATION  AND  EVALUATION 


We  first  applied  the  technique  described  earlier  bo  a set  of 
temperature- vapor  pressure  data  collected  over  a sample  of  solid  aluninun 
chloride  (3) . Hie  data  are  presented  in  Table  I. 

We  assumed  vapor  pressure  of  aluminum  chloride  could  be  represented 
by  an  equation  of  the  form 

log  P = | + B (34) 


where  the  variables  p and  T represent  pressure  in  Torr  and  temperature 
in  Kelvin,  respectively;  and  A and  B are  two  parameters  whose  values  are 
sought. 

We  solved  for  values  of  A and  B which  minimized  the  sum 


n 

I 

i=l  1 


(35) 


where  z*,  9z*/9A  and  Bz*/3B  may  be  found  from  Equations  (13)  and  (14). 
In  the  present  case 


Api  = Pi  - exp  [ (A/lb  + B)  In  (10)] 
AT±  = T±  - A/(log(Pi)  - B] 


(36) 


and 


9p'/9A  = p'  ln(10)/T 
3T'/9A  = l/[log(p)  - B] 

(37) 

9p'/9B  = p'  ln(10) 

9T'/9B  = A/[log(p)  - B]2 

The  distance  function,  z* , and  its  derivatives  were  calculated  for 
each  point  by  the  subroutine  beginning  at  step  6000  of  the  program  given 
in  the  next  section.  A careful  examination  and  evaluation  of  our 
particular  experimental  technique  and  apoaratus  led  to  the  selection  of 


V » • ■- 


*-  i 

I 
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I 
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TABLE  I.  VAPOR  PRESSURE  OF  SOLID  ALUMINUM  CHLORIDE 


Pressure,  Temperature, 
torr  °C 


3.77 

113.66 

M 3 

168.02 

158.28 

C 

3.80 

113.66 

C 

230.13 

163.40 

C 

3.92 

114.37 

M 

230.63 

163.40 

C 

5.52 

116 . 70 

C 

235.66 

163.40 

C 

5.53 

116.70 

M 

376.17 

169.90 

C 

6.30 

118.50 

M 

573.20 

176.80 

C 

6.32 

118.50 

C 

574.97 

176 . 80 

C 

7.14 

119.37 

M 

696.79 

179.66 

C 

7.25 

119 .37 

C 

716.93 

179.49 

C 

9.36 

122.87 

M 

808.52 

181.84 

C 

10.80 

125.61 

M 

934.15 

184.00 

C 

10.91 

124.23 

M 

943.10 

184.08 

C 

10 . 94 

124.23 

C 

952 . 31 

184.08 

C 

15.09 

127.97 

c 

1082.23 

186.25 

C 

16 . 26 

128.45 

M 

1157.09 

187.27 

C 

16.36 

128.45 

M 

1203.41 

187.87 

C 

16.47 

128.63 

C 

1344.23 

189.49 

C 

19.03 

130.92 

M 

1428.89 

190.49 

C 

33.81 

137.46 

C 

1461.56 

190.82 

C 

41.13 

140.42 

C 

1522.82 

191.44 

C 

58.15 

144 .53 

c 

1567.17 

191.95 

C 

88.42 

149.66 

c 

1641.14 

192.63 

C 

103.78 

152.49 

c 

1660.37 

192.79 

C 

106.05 

152.44 

c 

1705.21 

193.17 

C 

a M and  C refer  to  the  use  of  a McLeod  Gauge  or  a capacitame  manometer, 
respectively,  as  the  pressure  measuring  device. 
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6T  = 0. 1°K  and  6p  = 0.0002  p or  0.1  torr,  whichever  was  greater  for  each 
point. 

The  results  of  the  minimum  distance  least  squares  fit  are  shown  in 
Table  II.  The  fit  was  made  on  a Wang  2200  calculator  in  BASIC  language. 
The  program  is  listed  in  the  next  section.  Conventional  least  squares 
fits  also  were  made  of  these  data  for  comparison,  where  first  temperature 
then  pressure  was  made  to  be  the  independent  variable.  The  results  of 
all  three  least  squares  fits  are  summarized  in  Table  III. 

While  one  my  not  a priori  state  which  set  of  A and  B is  "best"  for 
these  data,  several  interesting  observations  can  be  mde.  The  value  of 
dP/dT  changes  from  0.34  to  107  over  the  range  of  the  data.  This  means 
that  a small  error  in  one  variable  propagates  a relatively  huge  error  in 
the  other  variable  at  one  end  of  the  data  set.  Therefore,  by  assuming 
either  one  of  the  variables  to  be  exact,  the  least  squares  fit  will 
reproduce  the  data  relatively  well  at  one  extreme  of  the  data  range  and 
relatively  poorly  at  the  other  extreme.  Furthermore,  it  my  be  seen  that 
the  data  points  are  on  the  average  only  0.18  torr  and  0.24K  from  the 
line  obtained  from  the  miniirarn  distance  fit.  If  temperature  were  assumed 
exact  (probably  the  conventional  assumption  in  least  squares  fitting  of 
data  to  Equation  34) , pressures  are  on  the  average  5.5  torr  in  error. 

In  order  to  examine  the  question  of  which  line  is  actually  "best," 
we  fit  a set  of  hypothetical  p,T  data  obtained  as  follows:  a "true" 

equation  was  assumed  for  which  A = -6000  and  B = 16.  Temperature  was 
stepped  at  equal  increments  from  104°C  to  198°C  to  yield  48  points,  and 
the  corresponding  pressures  calculated  from  the  equation.  These  "true" 
pressures  and  temperatures  were  then  randomized  byp=  p + 26p(0.5-R)  and 
T = T + (2)  (0.1)  (0.5-R) , where  6p  is  0.0002p  or  0.1,  whichever  is  greater, 
and  R is  a random  number  for  each  datum  such  that  0<R<1. 
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tart J*  xi.  MINIMUM  DISTANCE  LEAST  SQUARES  FIT  FOR  ALUMINUM  CHEORIDE 


VAPOR  PRESSURE  OF  SOLID  ALUMINUM  CHLORIDE 

-5891.065581615  +/-  12.16813381809 

15.8615816687  + /-  2 . 7 5 345 336 E-02 


ITN®  5 

SIGMA® 

2 . 989E+00 

DELTA  SIGMA® 

6. 689E-11 

NR  VARIABLE  NAME 

UNITS 

METRIC 

1 PRESSURE 

2 TEMP 

TORR 

K 

COMPUTED 
1 . 000E-01 

VARIABLE 

DELTA 

DELTA  PRIME 

ID  / Z* 

3.77 

-5 . 126E-01 

-4 . 521E-01 

M 

386.81 

1.401E+00 

1.654E-01 

4. 814E+00 

3.8 

-4.826  E-01 

-4 . 253E-01 

C 

386.81 

1. 314E+00 

1.561E-01 

4.530E+00 

3.92 

-6 . 468E-01 

-5 . 636E-01 

M 

387.52 

1. 683E+00 

2. 165E-01 

6 . 038E+00 

5.52 

-1 . 095E-01 

-8 . 7 7 8 E-02 

C 

389.85 

2. 200E-01 

4. 370E-02 

9. 806E-01 

5.53 

-9 . 9 54  E-02 

-7.9  74  E-0  2 

M 

389.85 

1.997  E-01 

3 . 9 7 3 E-02 

8. 909E-01 

6.3 

-3 . 05  7 E-01 

-2. 305E-01 

M 

391.65 

5 . 351E-01 

1. 317E-01 

2. 654E+00 

6.32 

-2 . 85  7 E-01 

-2.152  E-01 

C 

391.65 

4 . 994E-01 

1 . 2 31E-01 

2.480E+00 

7.14 

7. 165E-03 

5.137  E-0  3 

M 

392.52 

-1 . 140E-02 

-3.22  7 E-0  3 

6. 067E-02 

7.25 

1. 171E-01 

8.366  E-02 

C 

392.52 

-1 . 851E-01 

-5.29  4 E-0 2 

9 . 900E-01 

9.36 

-3. 206E-01 

-1.909  E-01 

M 

396.02 

3 . 890E-01 

1.57  3E-01 

2.474  E+00 

10 . 8 

-1 .449E+00 

-7 . 353E-01 

M 

398.76 

1.470E+00 

7.24  6 E-0 1 

1. 032E+01 

TABLE  II 


Continued 


10.91 

2.532E-02 

1. 349E-02 

M 

397 . 38 

-2 . 705E-02 

-1.263E-02 

1. 848E-01 

10.94 

5 . 532E-02 

2. 945E-02 

C 

397.38 

-5 . 902E-02 

-2 . 760E-02 

4 . 0 36E-01 

15.09 

1.262E-01 

4 . 84  8 E-0  2 

C 

401 . 12 

-9 . 969E-02 

-6. 140E-02 

7. 824E-01 

16.26 

6 . 790E-01 

2 . 4 35  E-0 1 

M 

401.6 

-5.078E-01 

-3.25  6 E-0 1 

4.067  E+00 

16.36 

7 . 790E-01 

2 . 7 84  E-0 1 

M 

401.6 

-5.809E-01 

-3. 733E-01 

4.65  7E+00 

16.47 

6.514E-01 

2.29  7E-01 

C 

401.78 

-4 . 808E-01 

-3 . 112E-01 

3. 868E+00 

19.03 

-1.241 E-0 1 

-3 . 5 2 9 E-0  2 

M 

404.07 

7 . 823E-02 

5.599  E-0  2 

6.618E-01 

33.81 

1. 115E+00 

1 . 3 70E-01 

C 

410.61 

-4 . 174E-01 

-3. 661E-01 

3 . 9 09  E+00 

41.13 

-2 . 851E-01 

-2 . 4 3 3 E-0  2 

C 

413.57 

8. 711E-02 

7.967E-02 

8. 331E-01 

58.15 

9 . 59  2 E-0 1 

4 . 549E-02 

C 

417 . 68 

-2 . 140E-01 

-2 . 0 38  E-0 1 

2 . 0 88  E+00 

88.42 

3 . 608E+00 

8. 188E-02 

C 

422.81 

-5 . 498E-01 

-5 . 373E-01 

5.435E+00 

103.78 

-1. 196E+00 

-1 . 926E-02 

C 

425.64 

1.530E-01 

1. 506E-01 

1.518E+00 

106.05 

1.465E+00 

2. 320E-02 

C 

425.59 

-1 . 859E-01 

-1.829  E-0 1 

1 . 844E+00 

168.02 

7 . 015E+00 

4.86  IE-02 

C 

431.43 

-5 . 860  E-0 1 

-5.819E-01 

5 . 8 39  E+00 

230.13 

-2.671E+00 

-9 . 800E-03 

C 

436.55 

1.621E-01 

1. 6 15E-01 

1 . 618E+00 

230.63 

-2. 171E+00 

-7.95  0 E-0  3 

C 

436.55 

1. 316E-01 

1. 311E-01 

1. 314E+00 

235.66 

2 . 858E+00 

1 . 0 2 5 E-0  2 

C 

436.55 

-1. 715E-01 

-1.70  8 E-0 1 

1. 711E+00 

376.17 

8.916E+00 

1. 351E-02 

C 

443.05 

-3 . 474E-01 

-3 . 4 68  E-0 1 

3.47 1E+00 

TABLE  II.  Continued 


573.2 

-1 . 4 11E+01 

-1 . 224  E-02 

C 

449.95 

3.627  E-01 

3.624E-01 

3. 625E+00 

574.97 

-1.234E+01 

-1.074  E-02 

C 

449.95 

3 . 168E-01 

3.165  E-01 

3 . 166E+00 

696.79 

-1.37 1E+01 

-1.2  26  E-02 

C 

452.81 

2 . 943E-01 

2.94 1 E-01 

2. 942E+00 

716.93 

1. 4 37E+01 

1. 339E-02 

C 

452.64 

-3.061E-01 

-3 . 058E-01 

3 . 060E+00 

808.52 

-1. 164E+01 

-1.06  7 E-02 

C 

454.99 

2. 180E-01 

2. 178E-01 

2 . 1 79  E+00 

934.15 

-1.007E+01 

-9.451E-03 

C 

457.15 

1. 652E-01 

1.651 E-01 

1.651 E+00 

943.1 

-6.043E+00 

-5 . 697E-03 

C 

457.23 

9.842  E-02 

9. 832E-02 

9 . 837E-01 

952.31 

3 . 166E+00 

3. 016E-03 

C 

457.23 

-5 . 134E-02 

-5. 129E-02 

5 . 1 31 E-01 

1082.23 

-9.690E+00 

-9 . 285E-03 

C 

459 . 4 

1. 386E-01 

1. 385E-01 

1. 385E+00 

1157.09 

-8.644E+00 

-8 . 370E-03 

C 

460.42 

1. 162E-01 

1. 161E-01 

1. 162E+00 

1203.41 

-7. 890E+00 

-7.687  E-0  3 

C 

461.02 

1.023E-01 

1.022E-01 

1.023E+00 

1344.23 

1.474E+00 

1. 468E-03 

C 

462.64 

-1. 731E-02 

-1. 729E-02 

1. 730E-01 

1428.89 

-1.522E+00 

-1 . 526E-03 

C 

463.64 

1. 688E-02 

1 . 6 86  E-0  2 

1 . 6 87  E-01 

1461.56 

1.069E+00 

1 . 0 7 7 E-0  3 

C 

463.97 

-1. 161E-02 

-1 . 160E-02 

1. 161E-01 

1522 . 82 

4.221E+00 

4. 28  3E-03 

C 

464 .59 

-4.41 7E-02 

-4 . 4 12  E-02 

4 . 4 14E-01 

1567.17 

-8. 346E-01 

-8.477E-04 

C 

465 . 1 

8. 49QE-03 

8 . 4 82  E-0  3 

8 . 4 86  E-0  2 

1641.14 

4 . 9 30E+00 

5 . 05  5 E-0  3 

C 

465 . 78 

-4 . 812E-02 

-4 . 807E-02 

4 . 809E-01 

1660.37 

7.715  E+00 

7 . 9 36  E-0  3 

C 

465.94 

-7 . 455E-02 

-7.44  7 E-02 

7. 451E-01 
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TABLE  II.  Continued 


1705.21  1.287  E + 01  1.333E-02  C 

466.32  -1.215E-01  -1.214E-01  1.215E+00 


NR  VARIABLE  NAME  RMS  DELTA 


RMS  DELTA  PRIME  UNITS 


1 PRESSURE 

2 TEMP 


5.963E+00  1.841E-01  TORR 

5.008  E-0 1 2.355E-01  K 


48  POINTS 


A and  B are  the  parameters  sought  in  Equation  (34) . The  values 

are  the  estimated  errors  in  A and  B,  and  are  calculated  from  the 
appropriate  diagonal  element  of  the  inverted  H matrix,  h|^,  using 
the  formula  = sigma  x^BS(|£l)  <i>  • 

Five  iterations  were  required.  At  each  iteration  the  value  of  the 
standard  deviation,  sigma  was  calculated.  Sigma  = (£(z*2)/tt) lj/2. 
Iterations  were  continued  until  the  change  in  sigma  was  smaller  than 
some  arbitrary  value.  The  convergence  criterion  in  the  present  case 

-7 

was  A sigma  < 1 x 10  . 

The  "metric"  in  pressure,  i.e. , 6p,  was  carputed  for  each  point  by 
the  subroutine  at  line  5000  in  the  program. 

The  pairs  of  numbers  are  on  the  first  line:  p,  Ap,  A'p,  device  code; 
on  the  second  line:  T,  AT,  A'T,  and  z*. 


Root  mean  square  values  are  given  for  Ap,  A'p,  AT,  and  A'T. 


TABLE  III.  COMPARISON  OF  LEAST  SQUARES  FITS  FOR 
SOLID  ALUMINUM  CHUDRIDE 


Minimum 
Distance  Fit 

Tenperature 

Exact 

Pressure 

Exact 

A 

-5891±12 

-5946+18 

-5923+14 

B 

15.9+2.8 

16.0+4.0 

15.913.2 

RMS  Ap,  torr 

6.0 

5.5 

6.7 

RMS  A'p,  torr 

0.18 

5.5 

- 

RMS  AT,  K 

0.50 

0.51 

0.47 

RMS  A'T,  K 

0.24 

— 

0.47 

The  same  three  types  of  least  squares  fits  were  made  on  this  set 
of  hypothetical  data,  and  the  results  are  svnnarized  in  Table  IV.  As 


before,  it  is  clear  that  the  pressures  and  temperatures  are  much  closer 
to  the  ndnii|m  distance  fit  line  than  to  either  of  the  other  lines.  It 
is  also  obvious  that  the  exact  ten^jerature  fit  reproduced  the  "true" 
equation  much  more  poorly  than  did  the  other  two  fits,  yet  this  is 
probably  just  the  kind  of  least  squares  fit  that  most  investigators 
would  choose  to  irake  if  they  were  to  select  one  independent  variable 
and  one  dependent  variable  from  Equation  (34) . The  apparent  agreement 
between  the  minimum  distance  fit  and  exact  pressure  fit  values  of  A and 
B is  not  surprising  for  the  randan  errors  in  pressure  propagate  errors  in 
tenperature  snaller  than  the  actual  6T  over  most  of  the  data  range. 
Another  way  of  expressing  this  is  to  say  that  pressure  behaves  more 
nearly  like  the  independent  variable  and  tenperature  the  dependent  one. 

The  fact  that  all  three  fits  of  hypothetical  data  shew  snaller 
RMS  A'x^  and  RMS  A'x^  than  do  the  corresponding  fits  of  real  data  merits 
some  oonments.  The  two  principal  reasons  are,  first,  Equation  (34) 
assumes  a constant  enthalpy  of  sublimation  for  aluminum  chloride  and 
ideal  gas  behavior  of  the  vapor.  Neither  assumption  is  entirely  valid. 
The  second  reason  is  that  the  random  errors  introduced  into  the  hypo- 
thetical data  were  constrained  to  not  exceed  <5p  and  6T,  respectively. 

The  errors  in  the  real  data  were  not  so  constrained  to  not  exceed  these 
"maximum  expected"  error  limits. 

These  tvro  sets  of  least  squares  fits  strongly  suggest  that  the 
minimum  distance  criterion  is  valid  and  desirable,  and  that  the  technique 
described  here  yields  results  at  least  as  good  as  do  conventional 
techniques,  and  may  do  considerably  better  than  conventional  techniques 
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TABLE  IV.  COMPARISON  OF  LEAST  SQUARES  FITS  FOR 
HYPOTHETICAL  p,T  DATA  SET 

-6000 

"True"  equation:  log  p = — ^ — + 16 


Mininun 
Distance  Fit 

Tanperature 

Exact 

Pressure 

Exact 

A 

-5999.5+2.3 

-5979. 2±5. 2 

-6002. 9±3. 1 

B 

16.00±5.3 

15.95+1.1 

16.01±8.8 

RMS  Ap,  torr 

2.0 

1.6 

2.4 

RMS  A'p,  torr 

0.02 

1.6 

- 

RMS  AT,  K 

0.12 

0.19 

0.12 

RMS  A'T,  K 

0.04 

- 

0.12 

if  an  urirfise  choice  of  which  variable  is  the  dependent  one  is  made  in 
the  latter  case.  When  one  considers  fitting  more  than  two  variables, 
the  minimum  distance  technique  becomes  even  more  attractive  as  it 

increasingly  difficult  to  single  out  a lone  dependent  variable 
for  conventional  fitting. 
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BASIC  LANGUAGE  IMPLEMENTATION  OF  FITTING  PROCEDURE 


The  mini"*™  distance  fitting  technique  described  earlier  has  been 
implemented  in  an  extended  BASIC  language  on  a Wang  2200S  calculator 
system.  The  calculator  had  the  general  I/O  and  advanced  programing 
options  and  8K  bytes  of  memory.  Output  was  produced  with  a model  2202 
plotting  typewriter. 


I 

[ 


- 

* - 

j 

r i 


■ 


Program  Concept 

The  program  was  designed  to  fit  functions  with  up  to  eight  variables 
and  one  to  eight  parameters.  The  program  performs  the  following  functions: 

1.  Determine  the  number  of  variables,  their  names,  the  units 
in  which  they  were  measured,  and  their  uncertainty  in  measurement. 

2.  Produce  on  request  a cassette  tape  containing  the  data  in 
a highly  packed  format  retaining  seven  decimal  digit  significance. 

3.  Ascertain  the  number  of  parameters  to  be  determined,  their 
nanes,  and  the  fractional  change  in  sigma  necessary  to  terminate  the 
fitting. 

4.  Generate  the  necessary  matrices  for  fitting  from  user 
written  routines  for  the  calculated  values  of  the  variables  fx'}  and 
partial  derivatives  }.  Reduction  to  unit  independent  quantities  is 
accomplished  using  the  uncertainties  in  measurement,  or  "metric,"  or, 
if  specified,  the  oonputed  uncertainty  in  measurement  of  each  variable 
at  each  point  as  supplied  by  a user  written  routine. 

5.  Invert  the  necessary  matrix  and  solve  for  improved  values 
of  the  parameters,  calculate  the  sigma  and  fractional  change  in  sigma. 
Thirteen  decimal  digits  significance  is  retained  during  these  calculations. 
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6.  Display  the  iteration  number,  number  of  points  used,  sigma, 
change  in  sigma,  and  the  adjusted  parameters  on  each  iteration.  Test 
the  fractional  change  in  sigma  to  determine  if  further  iterations  are 
necessary. 

7.  On  termination  of  fitting,  print  the  following  information: 

a.  Title  of  fit. 

b.  The  par  arte ters,  their  values  and  estimated  uncertainties 

c.  The  number  of  iterations,  final  value  of  sigma  and 

last  fractional  change  in  sigma. 

d.  Each  variable,  its  position  in  the  data  list,  the  units 

it  is  in,  and  its  metric. 

e.  On  request,  the  value  of  the  measured  variables,  their 
error  assuming  all  other  variables  correct  (Ax) , the  estimate  of  the  real 
error  in  the  variable  (A'x) , the  reduced  distance  of  the  point  from  the 
functional  surface  (z*) , and  a point  identification  of  up  to  six  alpha- 
numeric synbols. 

f.  The  root-mean-square  values  of  the  Ax  and  A'x  for 

each  variable  in  the  units  of  the  data. 

All  procedures  which  are  cannon  to  all  fitting  processes  have  been 
incorporated  into  the  main  program.  The  remaining  sections,  namely  the 
calculations  of  {x' } , {|^}  and  any  computed  metrics  must  be  supplied 
for  each  function  (or  set  of  data  for  the  computed  metrics)  by  the  user 

as  described  later. 

Program  Implementation 

The  program  used  for  fitting  the  examples  of  the  previous  section 
It  will  be  assumed  that  the  reader  is  familiar 


is  listed  in  Table  V. 


TABLE  V.  BASIC  PROGRAM  LISTING 


I 


I 

i 

! 


10  DIM  B(8) ,M(8),P1(8) ,S(8) ,T(8) ,X(8) , X0 ( 8 ) , H ( 8 , 8 ) , PO  (8 , 8 ) 

20  DIM  A$6,F$1,M$9,N1$1,N2$1,R$64,T$64, Y$1 

30  DIM  B$(4)64,C$(12)2,M$(8)1,T$(8)8,U$(8)8.X$(8)8 

40  SELECT  PRINT  005:  PRINT  HEX(03);:  INPUT  "TITLE  OF  FIT",T$ 

50  INPUT  "NUMBER  OF  VARIABLES" , N 

60  MAT  REDIM  X (N ) , X$ ( N ) 8 , U $ ( N ) 8 , XO  (N ) : BIN (Nl$ ) =N* 5 + 10 

70  PRINT  "NAMES  OF  VAR IABLE , UN  ITS , ME TRIC " : MAT  M=CON(N) 

80  FOR  1 = 1 TO  N:  INPUT  "[  ] " , X $ ( I ) ,U$  (I)  ,M(I) 

90  BIN(M$ (I) ) = I : IF  M(I)[=0  THEN  100:  BIN(M$(I))=0 

100  NEXT  I 

110  INPUT  "NEW  DATA  TAPE",Y$:IF  Y$[]"Y"  THEN  120:  GOSUB  1000 
120  SO , N0  = 0 : S=-9  9 : REWIND 

130  INPUT  "NUMBER  OF  PARAMETERS"  P : MAT  REDIM  T(P),T$(P)8 

140  PRINT  "ENTER  PARAMETER  NAME  S " : MAT  T=ZER:T=lE-7 :MAT  INPUT  T$ 
150  PRINT  "INPUT  DESIRED  DELTA  SIGMA  IF  OTHER  THAN  ";T;:INPUT  T 
160  FOR  1=1  TO  P : PRINT  T$(I);"  ESTIMATE" ;: INPUT  T(I):NEXT  I 
170  GOTO  300 

180  N0=N0+1 :N9 , S=0 : MAT  B=ZER(P):  MAT  H=ZER(P,P) 

190  GOSUB  'IOO(O) 

200  GOSUB  '100(1):  IF  N1]0  THEN  280 

210  BIN ( F$ ) =0  : GOSUB  ’200:  IF  VAL(F$)[]0  THEN  200 

220  GOSUB  '201:  Z1  = 0 : N9  = N9  + 1:  PRINT  HEX(01);N9:  MAT  P1=ZER(P) 

230  FOR  1=1  TO  N:  D=X ( I ) -XO (I ) : X= (M ( I ) /D) ! 2 : Z1=Z1+X 

240  FOR  K=1  TO  P:  PI (K) =P1 (K) -X*P0 (I , K) /D:  NEXT  K:  NEXT  I 

250  Z 1=  SQR ( 1 / Z 1 ) : S=S+Z1*Z1:  MAT  P1=(Z1!3)*P1 

260  FOR  1 = 1 TO  P:  B(  I)  =B  (I ) +Z  1 *P  1 (I ) : FOR  K=1  TO  P 

270  H(I,K)=H(I,K)+P1(I)*P1(K)  : NEXT  K:  NEXT  I:  GOTO  200 

280  FOR  1=1  TO  P:  FOR  K=1  TO  I:  H(  I , K)  =H  (K , I)  : NEXT  K:  NEXT  I 

290  MAT  P0= INV ( H) : MAT  S=P0*B:  MAT  T=T-S:  S=SQR(S/N9) 

300  S9=ABS (S-SO) / (S+S0)*2 

310  PRINT  : PRINT  HEX(03);"ITN  " ; NO ; " N=",N9,S,S9:  MAT  PRINT  T 

320  IF  ABS (S-SO) [T* (S+SO) /2  THEN  330:  SO=S:  GOTO  180 

330  STOP  '^CONVERGENCE  OCCURRED.  ENTER  continue  FOR  PRINTOUT" 

340  INPUT  "PRINT  ALL  DATA" , Y$ : SELECT  PRINT  213(80):  REWIND 

350  GOSUB  2000:STOP  "ENTER  continue  FOR  ADDN'L  OUTPUT":GOTO  340 

1000  N9  =0 : GOSUB  '100(2) 

1010  A$="END":  N9=N9+1 

1020  INPUT  "DATA  POINT:", A$:  IF  A$="END"  THEN  1060 

1030  MAT  INPUT  X:  Y$="Y":  INPUT  "DATA  OK",Y$ 

1040  IF  Y$  = "Y " THEN  1050:  PRINT  "RE-ENTER  ";A$;"  GOTO  1020 

1050  GOSUB  '100(3):  GOTO  1020 

1060  INPUT  "ALL  DATA  ENTERED", Y$:  IF  Y$[]"Y"  THEN  1010 

1070  GOSUB  '100(4):  PRINT  "DATAFILE  CLOSED":  RETURN 
1080  DEFFN ' 100 (0) : ON  Q GOTO  1120,1190,1200,1240 
1090  DATA  LOAD  ’’DATAFILE" 

1100  DATA  LOAD  BTB$  ( ) : IF  END  THLN  1180 

1110  Nl-0:  BIN(N2$)«1:  IF  Q=0  THLN  1170 

1120  IF  252-VAL (N2$) [VAL(N1$)  THLN  1100 
1130  $UNPACKB$ () [VAL(N2$) . VAL(N1$) -4 ] TO  R$ 

1140  IF  STR ( R $ , 1 , 3 ) [ ] "END"  THEN  1 1 5 0 : N 1 = 7 : BACKS PAC E IF:  RETURN 
1150  UNPACK(  + // . iHHHHHI  ! ! ! ! )STR(R$  , 7)  TO  X() 

1160  A$  = STR(R$ , 1,6)  : ADD(N2$,N1$) 

1170  RETURN 


! J 
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TABLE  V-  Continued 


1180 
1190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
2000 
2010 
2020 
2030 
2040% 
2050 
2060% 
2070% 
2080 
2090 
2100 
2110 
2120 
2130 
2140 
2150 
2160 
2170 
2180 
2190 
2200 
2210 
2220 
2230 
2240% 
2250% 
2260 
2270 
2280% 
2290% 
2300 
2310 
2320 
5100 
5110 
6000 
6010 
6020 
6030 
6040 
60  50 
6060 
6070 
6080 


PRINT 


PRINT 


GOSUB 


BACKSPACE  IF:  RETURN  CLEAR  : RETURN 

REWIND  : DATA  LOAD  "DATAFILE":  BIN(N2$)=1:  RETURN 

STR(R$, 1,6) =A$:  PACK (+#.######!!!!) STR (R$ , 7 ) FROM  X() 
$PACKB$ ( ) [VAL(N2$) ,VAL(N1$)  ] FROM  STR ( R$ , 1 , VAL (N 1 $ ) -4 ) 

ADD (N2$,N1$) : IF  2 5 2 - VAL ( N 2 $ ) [ VAL ( N 1 $ ) THEN  1230:  RETURN 

DATA  SAVE  BTB$ ( ) : BIN(N2$)=1:  RETURN 

A$  , STR  (R$  , 1,6)  ="END"  : PACK(+#  . //#//###  ! ! ! ! ) STR(R$  , 7 ) FROM  X() 

$PACKB$ ( ) [ VAL ( N2  $ ) , VAL(N1$) ] FROM  STR (R$ , 1 , VAL ( N 1 $ ) -4 ) 

DATA  SAVE  BTB$ ( ) : DATA  SAVE  END  : BACKSPACE  IF:  RETURN 

GOSUB  ' 100  (0) :N9  = 0 :MAT  B=ZER ( N)  : MAT  S=ZER(N):MAT  P1  = ZER(N) 
PRINT  HEX ( 0D0 AOAOA) ;T$ :HEX(0A) : FOR  1=1  TO  P 

PRINT  T$(I) ,T(I) ,"  +/-" ; S*SQR(ABS (PO (I , I) ) ) 

NEXT  I:  PRINT  : PRINTUSING  20 4 0 , NO  , S , S 9 : PRINT  : PRINT 

ITN  = ###  SIGMA = #.###!!!!  DELTA  SIGMA=  -#.###!!!! 
PRINTUSING  2060:  PRINT 
NR  VARIABLE  NAME  UNITS  METRIC 
//  ########  ########  ######### 

FOR  1=1  TO  N 

IF  VAL (M$ (I) ) =0  THEN  2100:  M$  = " COMPUTED " : GOTO  2110 
CONVERT  M ( I ) TO  M $,(//.###!!!!) 

PRINTUSING  2070 , I,X$ (I) ,U$ ( I) ,M$ : NEXT  I:  PRINT  : PRINT 
IF  Y$  [ 1 " Y " THEN  ’2130:  PRINTUSING  2240  : PRINT 
GOSUB  ’ 100 ( 1 ) : IF  N1[]0  THEN  2270:  BIN(F$)=0:  GOSUB  '200 

IF  VAL ( F $ ) [ ]0  THEN  2130 

Z1=0 : N9=N9+1:MAT  X0=X-X0:  FOR  1=1  TO  N 
Z 1=  Z 1 + (M ( I ) / XO ( I ) ) !2:  NEXT  I:  Z1  = SQR(1/Z1) 

FOR  1 = 1 TO  N:  B(I)  = (Z1*M(I) ) ! 2 / XO  ( I ) 

S(I)=S(I)+X0(I) ! 2 : P1(I)=P1(I)+B(I)!2 
IF  Y$[]"Y"  THEN  2260:  PRINT  X(I), 

IF  I [ ] 1 THEN  2210:  M$=A$ 

IF  I ] N THEN  2220:  CONVERT  Z1  TO  M$ ,(#.###!!! I ) 

IF  I 2 THEN  2230:  IF  I=N  THEN  2230:  M$="  " 

PRINTUSING  2250, XO (I) ,B(I) ,M$ 

VARIABLE  DELTA  DELTA  PRIME  ID/Z* 

-// . IHHf  !!!  ! -#.###  ! ! ! ! ########## 

NEXT  I:  IF  Y$ [ ] "Y"  THEN  2130:  PRINT  : GOTO  2130 

PRINT  HEX ( 0 A0A0A0A) : PRINTUSING  2280:  PRINT 

NR  VARIABLE  NAME  RMS  DELTA  RMS  DELTA  PRIME  UNITS 
# ########  # . ### ! ! ! ! #.###!!!!  ######## 
FOR  1=1  TO  N 

PRINTUSING  2290 , I,X$ (I) , SOR ( S ( I) /N9) , SQR(P1 (I) /N9) ,U$ (I) 
NEXT  I:  PRINT  : PRINT  N9  ; '*  POINTS":  RETURN 

M(J)=.0002*X(1) : IF  M(J)].l  THEN  5110:  M(J)=.l 

RETURN 
DEFFN' 200 

X(2)=X(2)+273. 15:  XO (1) =10 ! (T(2)+T(l)/X(2)) 

XO (2 ) =T( 1) / (LOG(X(l) ) /LOG(IO) -T(2)  ) 

FOR  J=1  TO  N : ON  VAL(M$(J))  GOSUB  5100:  NEXT  J:  RETURN 
BIN(F$)=7:  RETURN 

DEFFN ' 20 1 : PO ( 1 , 2 ) =X0 ( 1 ) *LOG ( 10 ) 

P0(1,1)=P0(1,2) / X ( 2 ) 

PO(2,1)=1/(LOG(X(1))/LOC(10)-T(2)) : P0(2,2)=T(1)*P0(2,1)!2 
RETURN 


UNITS 

######## 


I,X$(I)  ,S 


S0R(S(1)/N9),SQR(P1(I)/N9),U$(I) 
; " POINTS":  RETURN 


M ( J ) ] . 1 THEN  5110  : 


RETURN 


■T(1)*P0(2,1) ! 2 
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with  standard  Dartmouth  BASIC  and  the  extensions  provided  for  the  Wang 
system  used.  The  character  set  of  the  typewriter  used  required  three 
nonstandard  synbols  in  the  listing.  They  are 

Listing  Standard  Svirbol 


The  extensions  to  BASIC  used  may  be  found  in  the  System  2200S  Wang 
BASIC  Language  Programming  Manual  published  by  Wang  Laboratories,  Inc. , 
Tewksbury,  MA. 

Statements  10-30  allocate  space  for  matrices  and  alphanumeric 
variables.  The  symbols  have  been  kept  consistent  with  earlier  parts  of 
this  report.  Primed  quantities  have  had  a zero  appended  and  starred 
c3uantities  have  had  a one  appended  since  primes  and  asterisks  are  not 
legal  parts  of  a BASIC  variable  (e.g. , x^  becomes  Xtf(I)  and  z*  becomes 
Zl). 

Statements  40-170  perform  the  functions  indicated  in  1,  2 and  3 
above.  The  sequence  of  statements  to  generate  a tape  appears  at  1000- 


1070. 


Statements  180-330  perform  the  functions  described  in  4,  5,  and  6 
above.  These  operations  cure  repeated  using  the  data  from  cassette  until 
convergence  of  the  values  is  achieved. 

Statements  340—350  with  the  associated  subroutine  sequence  in 
statements  2000-2320  perform  the  output  as  specified  in  7 above. 

Statements  1080-1260  are  a sequence  of  statements  for  cassette 
handling — to  read  and  write  the  data  on  tape.  The  argument  Q passed  to 
the  subroutine  determines  the  operation  to  be  performed. 


Value  of  Q (integer) 

Q<1  or  Q>4  Rewind  tape  and  open  "DATAFILE"  for  reading 

Q = 1 Read  a data  point  from  tape 

Q = 2 Rewind  tape  and  open  "DATAFILE"  for  writing 

Q = 3 Write  a data  point  to 

Q = 4 Write  any  data  points  stored  in 

buffer  to  tape  and  close  "DATAFILE"  on  tgpe 
The  following  three  sets  of  statements  are  user-supplied  and  tailored 


Q = 4 


to  each  fit: 

Statements  5100-5110  - compute  a metric  for  the  first  variable 
Statements  6000-6040  - ocrpute  the  values  of  the  primed 
variables  {x* } 

Statements  6050-6080  - conpute  the  partial  derivatives  } 

Program  Modification  for  Fitting  Other  Functions 

To  fit  another  function  the  statements  from  6000  on  must  be  replaced. 

A new  subroutine  to  calculate  {xV } must  be  written  (DEFEN  ’200) , and  a 
new  subroutine  to  calculate  the  necessary  partied  derivatives,  } 

must  replace  the  old  DEFFN  '201. 

In  addition,  if  the  method  of  conputing  a metric  for  any  variable 
is  changed,  suitable  changes  in  statement  6030  and  the  subroutine  at  5100 
would  have  to  be  made  also. 

1.  Calculation  of  the  {x^}  is  done  in  subroutine  DEFEN  '200. 

When  the  subroutine  DEFFN  *200  is  entered  the  measured  values  for  a 
point  cue  available  in  X(I) , where  I = 1 to  N (N  is  the  number  of  variables) , 
and  the  current  parameter  estimates  are  available  in  T(J)  where  J = 1 to 
P (P  is  the  number  of  parameters) . The  values  of  the  {xV } are  to  be 


calculated  and  stored  in  X0(I) , I = 1 to  N.  In  the  listing  the 
pressure  p'  = X0(1)  is  calculated  using  the  parameters  T(l) , T(2)  and  the 
temperature  t = X(2) . The  temperature  t ' = X0(2)  is  calculated  from  the 
parameters  T(l) , T(2)  and  the  measured  pressure  p = x(2) . 

2.  Calculation  of  partial  derivatives  } is  done  in 

subroutine  DEFETJ  '201.  The  partial  derivative  of  every  variable  with 
respect  to  each  parameter  must  be  calculated  and  stored  in  an  N x P array. 
The  partial  derivative  (|~1)  is  calculated  and  stored  in  P0(l,2) , and 


<!?>  “ W(2,l) 


3.  Calculation  of  computed  metrics  is  done  by  a normal  sub- 
routine internally  in  the  marked  subroutine  DEFFN  '200.  The  matrix 
element  M$(l)  contains  a zero  for  fixed  metric  variables,  or  the  index 
of  the  variable  for  variables  whose  metric  is  to  be  computed.  This 


index  may  be  used  to  call  the  correct  subroutine  as  in  the  example  of 
statement  6030  and  the  subroutine  that  statement  calls  at  statements 
5100-5110.  The  computer  metric  for  variable  I is  stored  in  M(I) . 


Operation 


Write  the  necessary  subroutines  and  metric  calculating  procedures. 
Insure  that  the  statement  numbers  for  these  routines  begin  with  3000 
or  a larger  number  (9999  is  the  largest  statement  number  allowed  on  the 
Wang  2200) . Save  the  program  cm  tape,  if  desired. 


(a)  Mount  the  tape  containing  QENFIT  and  rewind  it. 

(b)  Enter  LOAD  "GENFIT"  i where  i stands  for  RETURN  (EXEC) . 

(c)  Enter  RUN  i. 

C 

(D)  The  program  will  respond  with:  TITI£  OF  FIT? 

(d)  Enter  any  character  string  up  to  64  characters  long  which 
does  not  contain  a conma,  then  i. 
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(F)  The  program  will  respond  with:  NAMES  CF  VARIABLE,  UNITS, 

METRIC,  <>? 

(f)  For  each  variable,  when  the  <>?  is  displayed,  enter  the 
requested  information,  separated  by  commas,  then  4-.  After  the  last 
entry  program  response  is: 

(G)  NEW  DATA  TAPE? 

(g)  Enter  any  word  starting  with  a Y and  ending  with  a 4 

for  an  affirmative  answer.  Any  other  response  will  be  interpreted  as  NO  4. 
If  a data  tape  is  already  in  existence,  mount  and  rewind  it,  proceed  to 
(M)  belcw.  If  a new  tape  is  to  be  created,  insure  a tape  with  a data 
file  header  is  mounted  and  rewound.  (Note:  A data  file  header  may  be 
created  on  a new  tape  by  mounting  and  rewinding  the  tape  with  the  write 
tabs  on,  then  entering  the  oortmand  DATASAVE  OPEN  "DATAFILE"  4.) 

(H)  Program  response  after  affirmative  response  to  new  tape 
is  DATA  POINT:? 

(h)  Enter  any  point  identifier  (up  to  6 letters)  then  4 or 
if  last  data  point  has  been  entered,  enter  END  4. 

(I)  For  program  response  to  END  4 see  (K)  below.  For  any 
other  identifier:  ? 

(i)  Enter  variable  values  in  order  separated  by  ccnmas  or 
carriage  returns. 

(J)  Program  response  after  last  value:  DATA  OK? 

(j)  Enter  4 or  Y4  if  data  is  correct.  If  not  all  right, 
enter  any  string  of  characters  not  starting  with  a Y.  Program  response 
if  data  is  OK  is  to  return  to  output  (H) . If  data  is  not  OK,  data  will 
not  be  put  on  tape  and  return  to  (H)  for  re-entry  of  data  point. 


(K)  Program  response:  ALL  DATA  ENTERED? 

(k)  If  all  entered,  enter  Y4-  or  + ; more  data  to  be  entered 
any  response  not  beginning  with  Y will  return  program  to  (H)  above. 

(L)  Program  response  to  an  affirmative  to  all  data  entered 
is  to  close  the  data  file  and  print,  DATAFILE  CLOSED;  no  response  is 
required. 

(M)  Program  will  ask:  NUMBER  OF  PARAMETERS? 

(m)  Enter  the  nunber  of  parameters  the  program  is  designed 
to  fit,  (this  nurrber  is  called  P) , then  +. 

(N)  Program  will  respond:  ENTER  PARAMETER  NATES  ? 

(n)  Enter  the  parameter  names  in  order  separated  by  comas, 

then  I. 

(O)  Program  will  respond:  INPUT  DESIRED  DELTA  SIGMA  IF 

OTHER  THAN  1 . 00000000E-07? 

(o)  If  a different  value  of  the  fractional  change  in  sigma 
to  terminate  the  fit  is  desired  enter  it,  then  1,  otherwise  just  1. 

(P)  Program  will  respond  with  requests  for  each  of  the 

parameters ESTIMATE?  (Where is  the  parameter  name.) 

(p)  Enter  the  estimate  for  each  parameter  as  requested, 
followed  by 

(Q)  After  last  estimate  is  entered,  the  program  will  proceed 

to  iteratively  fit  the  data  on  the  casette  until  the  fractional  change 
in  sigma  is  smaller  than  the  entered  delta  sigma.  Initially  and  after 
each  iteration,  the  iteration  nutber,  nurrber  of  points  used,  sigma, 
the  fractional  change  in  sigma,  and  all  of  the  current  parameters  are 
displayed.  On  convergence  (fractional  change  in  sigma  < delta  signa) 
the  program  will  respond:  STOP:  ODNVERfiENCE  OCCURRED.  ENTER  CONTINUE 


i 


(q)  Enter  CONTINUE  + 

(R)  Program  will  respond  PRINT  ALL  DATA? 

(r)  If  all  data  is  to  be  printed  with  error  indicators 
respond  Y4-;  if  not  and  only  a summary  is  desired  give  any  response  not 
starting  with  Y and  only  the  sunrary  will  be  printed. 

(S)  Program  terminate  with  the  message  STOP  ENTER  CONTINUE 
FOR  ADDN'L  OUTPUT 

(s)  If  further  output  is  desired,  enter  CONTINUE  + and 
program  will  resune  at  (R)  above. 

Index  of  Variables  Used  in  Program 
Single  variables  (numeric) 

D (temporary)  Used  to  store  Ax^  = (x^-x^ ' ) 

I (temporary)  Used  as  an  index  in  FOR  loops 

J (temporary)  Used  as  an  index  in  FOR  loops 

K (temporary)  Used  as  an  index  in  FOR  loops 

N - Number  of  variables 
N0  - Iteration  number 

Nl  - End  of  file  indicator  for  datafile  = 0:no  EOF;  jt  0:EQF 
N9  - Number  of  points  entered  into  fit  (incremental  as  each  point 
is  processed) 

P - Number  of  parameters 

2 

S - Sigma;  used  to  collect  E(z*)  during  each  iteration 
SO  - Sigma  of  previous  iteration 

S9  - Fractional  change  in  sigma  on  last  iteration  ABS (S-SO) *2 

?S+S0l 

T - Delta  sigma;  termination  of  fit  indicated  by  S9<T 
X (temporary)  Used  to  calculate  z* 
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Zl  - z*;  the  distance  from  the  functional  surface  in  the  metric  spare 

Single  variables  (alphanumeric) [final  digits  indicate  length  in  characters] 
A$6  (temporary)  Point  identification  for  current  point 
F$1  (temporary)  Flag  used  to  indicate  valid  point  F$=0,  valid  point; 

F$^0,  invalid  point  program  will  not  use  for  fitting 
M$9  (tenporary)  Used  to  display 

Nl$l  - String  length  indicator  for  packing  and  unpacking  data  points 

N2$l  - Pointer  for  packing  and  unpacking  data  points 

R$64  (tenporary)  String  variable  used  to  pack  and  unpack  data  points 

T$64  - String  variable  used  to  hold  title  of  fit 

Y$1  - Response  indicator  for  various  portions  of  program  (YES/N3) 

Array  variables  (numeric) 

B(  ) - Used  for  error  vector  during  fitting;  used  to  store  (A'x) 
during  print  out 

M(  ) - Used  to  store  metrics  for  variables 

r3z*i  2 

Pl(  ) - Used  to  store  { gg— } during  fitting;  E(A'x)  during  printout 

1o 

S(  ) - Used  to  store  E (Ax)  during  printout 

T(  ) - Used  to  store  the  parameter  values 

X(  ) - Used  to  store  the  measured  values  for  each  point 

X0(  ) - Used  to  store  calculated  values  for  each  point  during  fit; 

used  to  store  {Ax}  during  printout 

H(  , ) - Curvature  matrix  for  least  squares  fit 

PO  ( , ) - Used  to  store  partial  derivatives  {-3^}  during  initial 

k 

part  of  fit;  used  to  store  the  error  matrix  (inverse  of  H) 


during  later  part  of  fit 


h 
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Array  variables  (alphanumeric) [final  digits  are  lengths  of  each  element 
in  characters] 

B$(  )64  - Buffer  used  to  store  data  for  tape  operations 

M$(  )1  - Used  to  store  ccirputed  metric  indices 

T$(  ) 8 - Used  to  store  names  of  parameters 

U$(  )8  - Used  to  store  name  of  units  for  variables 

X$(  ) 8 - Used  to  store  name  of  variables 


' l 
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